#!/usr/bin/perl

use strict;
use PGPLOT;

# get a list of all the spectra that we have
my $working_directory="/data/MENTOK_2/ste616/flux_monitoring/staging_area";
my $list_command="find ".$working_directory." -name '*.spectrum'";
print "EE executing $list_command\n";
open(LIST,"-|")||exec $list_command;
my @spectrum_files;
while(<LIST>){
    chomp;
    push @spectrum_files,$_;
}
close(LIST);

# get a list of shortest baseline lengths as a function of date
my @epoch_date;
my @epoch_short_baseline;
open(LOG,"uranus_1934_replot.run");
my $current_date;
while(<LOG>){
    chomp;
    my $line=$_;
    if ($line=~/^\+\+\sfound\sUranus/){
	my $date=$line;
	$date=~s/^(.*)time\:\s(.*?)\s(.*)/$2/;
	$current_date=$date;
    } elsif ($line=~/^DD\sshortest\sbaseline/){
	my $shortbase=$line;
	$shortbase=~s/^(.*)\=\s(.*?)\s(.*)/$2/;
	push @epoch_date,$current_date;
	my $shortbase_string=sprintf "%.0f",$shortbase;
	push @epoch_short_baseline,$shortbase_string;
    }
}
close(LOG);

my $dolog=1;

# get the MIRIAD model fluxes for 1934-638
open(OLDMODEL,"oldmodel_1934_fluxes.out");
my @oldmodel_1934_frequency;
my @oldmodel_1934_flux;
while(<OLDMODEL>){
    chomp;
    my @els=split(/\s+/);
    if ($dolog==1){
	push @oldmodel_1934_frequency,log($els[0])/log(10.0);
	push @oldmodel_1934_flux,log($els[1])/log(10.0);
    } else {
	push @oldmodel_1934_frequency,$els[0];
	push @oldmodel_1934_flux,$els[1];
    }
}
close(OLDMODEL);

# go through these files and make some plots
my @j1934_frequency;
my @j1934_flux;
my @j1934_date;
my @j1934_colour;
my @j1934_symbol;
my @uranus_frequency;
my @uranus_flux;
my @uranus_date;
my @uranus_colour;
my @uranus_symbol;
open(OUT,">combined_spectra.txt");
my $outfilen=-1;
for (my $i=0;$i<=$#spectrum_files;$i++){
    $spectrum_files[$i]=~/^\/(.*)\/(.*?)\/([^\/]*)$/;
    my $date_frequency=$2;
    my $source_spectrum=$3;
    $date_frequency=~/^(.*?)\_(.*)/;
    my $this_date=$1;
    $source_spectrum=~/^(.*?)\.(.*)/;
    my $this_source=$1;
    if ($this_source eq "1934-638"){
	$outfilen++;
    }
    open(SPECTRUM,$spectrum_files[$i]);
    while(<SPECTRUM>){
	chomp;
	my $line=$_;
	my @els=split(/\s+/,$line);
	if ($this_source eq "uranus"){
	    if ($dolog==1){
		push @uranus_frequency,log($els[1])/log(10.0);
		push @uranus_flux,log($els[2])/log(10.0);
	    } else {
		push @uranus_frequency,$els[1];
		push @uranus_flux,$els[2];
	    }
	    push @uranus_date,$this_date;
	} elsif ($this_source eq "1934-638"){
	    print OUT "$outfilen $line\n";
	    if ($dolog==1){
		push @j1934_frequency,log($els[1])/log(10.0);
		push @j1934_flux,log($els[2])/log(10.0);
	    } else {
		push @j1934_frequency,$els[1];
		push @j1934_flux,$els[2];
	    }
	    push @j1934_date,$this_date;
	}
    }
    close(SPECTRUM);
}
close(OUT);

my @done_dates;
my $current_colour=0;
my $current_symbol=2;
for (my $i=0;$i<=$#uranus_date;$i++){
    # check if this date has been done before
    my $isdone=0;
    for (my $j=0;$j<=$#done_dates;$j++){
	if ($done_dates[$j] eq $uranus_date[$i]){
	    $isdone=1;
	    last;
	}
    }
    if ($isdone==1){
	next;
    }
    $current_colour++;
    if ($current_colour==15){
	$current_colour=1;
	$current_symbol++;
    }
    for (my $j=0;$j<=$#uranus_date;$j++){
	if ($uranus_date[$j] eq $uranus_date[$i]){
	    $uranus_colour[$j]=$current_colour;
	    $uranus_symbol[$j]=$current_symbol;
	}
    }
    for (my $j=0;$j<=$#j1934_date;$j++){
	if ($j1934_date[$j] eq $uranus_date[$i]){
	    $j1934_colour[$j]=$current_colour;
	    $j1934_symbol[$j]=$current_symbol;
	}
    }
    push @done_dates,$uranus_date[$i];
}

# find the maximum and minimum values for all the axes
my ($minfreq,$maxfreq)=&minmax(@uranus_frequency,@j1934_frequency);
my ($j1934_minflux,$j1934_maxflux)=&minmax(@j1934_flux);
my ($uranus_minflux,$uranus_maxflux)=&minmax(@uranus_flux);

print "II 1934 maximum flux = ".$j1934_maxflux." @ ";
for (my $i=0;$i<=$#j1934_flux;$i++){
    if ($j1934_flux[$i]==$j1934_maxflux){
	print $j1934_date[$i]." @ ".$j1934_frequency[$i]."\n";
	last;
    }
}

print "II Uranus maximum flux = ".$uranus_maxflux." @ ";
for (my $i=0;$i<=$#uranus_flux;$i++){
    if ($uranus_flux[$i]==$uranus_maxflux){
	print $uranus_date[$i]." @ ".$uranus_frequency[$i]."\n";
	last;
    }
}

# it's plotting time
#pgopen("41/xs");
pgopen("uranus_1934_flux_calibration_gooddata_fixedmodel4_uranusshortbase.ps/cps");
my ($pgvp_xmin,$pgvp_xmax,$pgvp_ymin,$pgvp_ymax);
pgqvp(0,$pgvp_xmin,$pgvp_xmax,$pgvp_ymin,$pgvp_ymax);
#print "DD PGPLOT viewport limits: x = {".$pgvp_xmin.",".$pgvp_xmax."} ".
#    "y = {".$pgvp_ymin.",".$pgvp_ymax."}\n";
# write the axis labels now
pgswin(0,1,0,1);
pglab("Frequency (GHz)","Flux (Jy)","");
my ($vp_ymin,$vp_ymax,$vp_yheight);
my ($vp_plot_xmin,$vp_plot_xmax,$vp_text_xmin,$vp_text_xmax);
$vp_yheight=($pgvp_ymax-$pgvp_ymin)/2.0;
$vp_ymax=$pgvp_ymax;
$vp_ymin=$vp_ymax-$vp_yheight;
$vp_plot_xmin=$pgvp_xmin;
$vp_text_xmax=$pgvp_xmax;
$vp_plot_xmax=$pgvp_xmax-($pgvp_xmax-$pgvp_xmin)/7.0;
$vp_text_xmin=$vp_plot_xmax;
# 1934 first, at the top of the plot
pgsvp($vp_plot_xmin,$vp_plot_xmax,$vp_ymin,$vp_ymax);
pgswin($minfreq,$maxfreq,$j1934_minflux,$j1934_maxflux);
if ($dolog==1){
    pgbox("BCLTS",0.0,0,"BCLNTS",0.0,0);
} else {
    pgbox("BCTS",0.0,0,"BCNTS",0.0,0);
}
my $printdate;
my $maxcolour=1000;
for (my $i=0;$i<=$#j1934_frequency;$i++){
    if ($j1934_colour[$i]<$maxcolour){
	if ($j1934_date[$i] ne $printdate){
	    print "II colour ".$j1934_colour[$i]." is from ".$j1934_date[$i]."\n";
	    $printdate=$j1934_date[$i];
	}
	pgsci($j1934_colour[$i]);
	pgpt1($j1934_frequency[$i],$j1934_flux[$i],$j1934_symbol[$i]);
    }
}
# plot the MIRIAD 1934 model as a line
pgsci(2);
pgline($#oldmodel_1934_frequency+1,\@oldmodel_1934_frequency,\@oldmodel_1934_flux);
# plot some more lines from RJS 12mm 1934 memo
# the conventional fit
my @conventional_fit;
my @sault_fit;
my @new_fit;
for (my $i=0;$i<=$#oldmodel_1934_frequency;$i++){
    if ($dolog==1){
	$conventional_fit[$i]=-30.7667+26.4908*($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))-
	    7.0977*(($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))**2)+
	    0.605334*(($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))**3);
	$sault_fit[$i]=-202.6259+149.7321*($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))-
	    36.4943*(($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))**2)+
	    2.9372*(($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))**3);
	$new_fit[$i]=26.7035-15.6502*($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))+
	    3.23542*(($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))**2)-
	    0.242125*(($oldmodel_1934_frequency[$i]+log(1000.0)/log(10.0))**3);
    } else {
	$conventional_fit[$i]=10**(-30.667+26.4908*log($oldmodel_1934_frequency[$i])/log(10.0)-
				   7.0977*((log($oldmodel_1934_frequency[$i])/log(10.0))**2)+
				   0.605334*((log($oldmodel_1934_frequency[$i]/log(10.0)))**3));
    }
}
pgsci(3);
pgline($#oldmodel_1934_frequency+1,\@oldmodel_1934_frequency,\@sault_fit);
pgsci(4);
pgline($#oldmodel_1934_frequency+1,\@oldmodel_1934_frequency,\@conventional_fit);
pgsci(6);
pgline($#oldmodel_1934_frequency+1,\@oldmodel_1934_frequency,\@new_fit);
# label the 1934 plot
my $label_1934_x=($maxfreq-$minfreq)/20.0+$minfreq;
my $label_1934_y=($j1934_maxflux-$j1934_minflux)/10.0+$j1934_minflux;
pgsci(1);
pgtext($label_1934_x,$label_1934_y,"1934-638");
# Uranus now, at the bottom of the plot
pgsci(1);
$vp_ymax-=$vp_yheight;
$vp_ymin-=$vp_yheight;
pgsvp($vp_plot_xmin,$vp_plot_xmax,$vp_ymin,$vp_ymax);
pgswin($minfreq,$maxfreq,$uranus_minflux,$uranus_maxflux);
if ($dolog==1){
    pgbox("BCLNTS",0.0,0,"BCLNTS",0.0,0);
} else {
    pgbox("BCNTS",0.0,0,"BCNTS",0.0,0);
}
for (my $i=0;$i<=$#uranus_frequency;$i++){
    if ($uranus_colour[$i]<$maxcolour){
	pgsci($uranus_colour[$i]);
	pgpt1($uranus_frequency[$i],$uranus_flux[$i],$uranus_symbol[$i]);
    }
}
# label the Uranus plot
my $label_uranus_x=($maxfreq-$minfreq)/20.0+$minfreq;
my $label_uranus_y=$uranus_maxflux-($uranus_maxflux-$uranus_minflux)/10.0;
pgsci(1);
pgtext($label_uranus_x,$label_uranus_y,"Uranus (fixed model)");
# write the text dates
pgsci(1);
pgsvp($vp_text_xmin,$vp_text_xmax,$pgvp_ymin,$pgvp_ymax);
pgswin(0,100,0,100);
my @done_colour;
my @done_symbol;
my $symbol_xloc=10;
my $text_xloc=18;
my $symbol_yloc=95;
my $symbol_ydiff=3;
pgsch(0.85);
for (my $i=0;$i<=$#j1934_colour;$i++){
    my $alreadydone=0;
    for (my $j=0;$j<=$#done_colour;$j++){
	if (($j1934_colour[$i]==$done_colour[$j])&&
	    ($j1934_symbol[$i]==$done_symbol[$j])){
	    $alreadydone=1;
	    last;
	}
    }
    if ($alreadydone==0){
	pgsci($j1934_colour[$i]);
	pgpt1($symbol_xloc,$symbol_yloc,$j1934_symbol[$i]);
	my $this_length;
	for (my $j=0;$j<=$#epoch_short_baseline;$j++){
	    if ($epoch_date[$j] eq $j1934_date[$i]){
		$this_length=$epoch_short_baseline[$j];
		last;
	    }
	}
	pgtext($text_xloc,$symbol_yloc-1,$j1934_date[$i]." (".$this_length."m)");
	$symbol_yloc-=$symbol_ydiff;
	push @done_colour,$j1934_colour[$i];
	push @done_symbol,$j1934_symbol[$i];
    }
}
pgclos();

sub minmax {
    # get the minimum and maximum values in an array
    my @array=@_;

    my $minval=$array[0];
    my $maxval=$array[0];
    for (my $i=0;$i<=$#array;$i++){
	if ($array[$i]<$minval){
	    $minval=$array[$i];
	}
	if ($array[$i]>$maxval){
	    $maxval=$array[$i];
	}
    }
    return ($minval,$maxval);
}

